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We study an 517(2) symmetric spin-3/2 model on the z = 3 Bethe lattice using the infinite 
Time Evolving Block Decimation (iTEBD) method. This model is shown to exhibit a rich phase 
diagram. We compute the expectation values of several order parameters which allow us to identify 
a ferromagnetic, a ferrimagnetic, a anti-ferromagnetic as well as a dimerized phase. We calculate 
the entanglement spectra from which we conclude the existence of a symmetry protected topological 
phase that is characterized by S — 1/2 edge spins. Details of the iTEBD algorithm used for the 
simulations are included. 



I. INTRODUCTION 



Quantum antiferromagnets are a challenging and inter- 
esting topic in condensed matter physics. The interplay 
between strong correlations, quantum fluctuations, and 
frustration can yield very exciting new phases with often 
unexpected properties. One dimensional quantum spin 
chains have been proven to be very useful to understand 
many of these interesting phases. For these systems, very 
powerful analytical and numerical methods exist. A ma- 
jor advance on the numerical side was the introduction 
of the Density-Matrix-Renormalization Group (DMRG) 
method which allows for the efficient simulation of one- 
dimensional systems^ One of the first successes of the 
then new DMRG method was its use to prove the fa- 
mous Haldane conjecture for integer spin chains numer- 
ically. The Haldane conjecture states that the Heisen- 
berg antiferromagnetic (HAF) chain with integer spin S 
has a nonzero excitation gap and exponentially decaying 
spin correlation functions (while spin chains with half- 
integer spin are gapless)!^ For odd integer spins, the 
Haldane phase is an example of a so-called symmetry 
protected topological phase (SPTP) which is character- 
ized by S = 1/2 edge spins. This kind of phase cannot 
be characterized by symmetry breaking but instead by 
using cohomology theory- ^ 

In this paper we make use of an extension of the DMRG 
algorithm to the Bethe lattice and study the phase dia- 
gram of a general SU{2) symmetric S = 3/2 spin model 
on a Bethe lattice with coordination number z = 3. This 
model has a special point in this phase diagram, the 
so-called Affleck-Kennedy-Lieb-Tasaki (AKLT) point, at 
which the ground state is known exactlypB The AKLT 
wave function for this model is a quantum paramagnet 
with exponentially decaying correlation functions 
We argue that the AKLT point extends to a phase which 
is similar to the Haldane phase in one-dimensional spin 
chains in that it is characterized by a fractionalized edge 
spin (the precise meaning of an edge spin in this con- 
text is explained Sec. IV). Beside the Haldane phase, the 



well as a dimerized phase. Even though there exist ex- 
perimental systems, such as dendrimers that realize the 
tree structure^ related to the Bethe lattice considered 
in this work, we are mainly interested in this system be- 
cause of its theoretical nature. The 5 = 3/2 model has 
already been shown to exhibit a very complex phase di- 
agram on the mean field level and in one-dimensional 
systems!^^ The main goal is to present a conclusive 
phase diagram of the model utilizing recently introduced 
algorithms which allow an efficient simulation of quan- 
tum spin systems on Bethe lattices. We use a descen- 
dant of Vidal's infinite time-evolving block decimation 
algorithm (iTEBD) 18 adapted to the tree like structure. 
The iTEBD method as well as the DMRG have already 
successfully been applied to reproduce, e.g., the phase 
diagram of the transverse field Ising model and the spin 
1/2 XXZ model on the z = 3 Bethe latticePHm We in- 
clude details of the algorithm in this paper and discuss a 
number of improvement that make it more stable. 

This paper is organized as follows: in Sec. [H] we in- 
troduce the model and discuss some of its basic prop- 
erties leading up to section Sec. |III| where we describe 
the method we used to obtain the results. In Sec. IIVI we 
first take a closer look at symmetry protected topological 
phases before we present the results in Sec. [V] The key 
points of this paper are summarized again in Sec. |VI| 



II. MODEL 

Throughout this paper we consider the following 
nearest-neighbor spin-3/2 model Hamiltonian 

H = Y^aS t - S l+1 + ■ S t+1 ) 2 + 7 (Si ■ § i+1 ) 3 (1) 



where a different parametrization of the Hamiltonian is 
given by 



model is shown to exhibit different magnetic phases as 
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with <p £ [— n, 7r] and € [—§,§]• The symmetries of 
this model include translation, spatial inversion, SU(2), 
and time reversal (TR). This model is known to ex- 
hibit an AKLT-likc wavefunction at the point a = 1.0, 
P = Ml' T = !H3 when placed on a z = 3 Bethe lat- 
tice. The AKLT state has symmetry protected S — 1/2 
edge spins which are discussed in detail in Sec. |IV| Fur- 
thermore this model exhibits a SU(A) symmetry at four 
points in the phase diagram conn ected by a SO(5) sym- 
metric line which is given b ^ 15 * 16 * a = — ^(31Jo + 23J2), 
P = i(5J + 17J 2 ), and 7 = ^(J + J 2 ) with J , J 2 > 0. 
We place this Hamiltonian on the Bethe lattice with co- 
ordination number z — 3. The Bethe lattice and its finite 
counterpart, the Cay ley tree, has first been used in sta- 
tistical mechanics PSES More recently it has also proved 
to be a highly instructive testing ground for tensor net- 
work methods as it is loop-free thus removing one of the 
major sources of entanglement. Additionally, the Bethe 
lattice is self-similar, enabling the application of efficient 
infinite-system methods. This lattice is infinite by defini- 
tion and thus there are no surface effects.^ Note that the 
thermodynamic limit of the Cayley tree and the Bethe 
lattice are not equivalent. This inequivalence is rooted in 
the large number of surface sites contained in any Cayley 
treeP^ Whereas in most systems the ratio of boundary 
to bulk sites reduces to zero for large systems, it remains 
finite for the Cayley tree. The finite tree is thus domi- 
nated by the boundary conditions, making it unsuitable 
to study the model's properties in the thermodynamic 
limit. 



III. TENSOR PRODUCT STATE BASED 
SIMULATIONS ON A TREE LATTICE 



is the "physical" index, e.g., enumerating the spin states 
on each site, and A is a x x X; real, diagonal matrix. 
Ground states of one dimensional gapped systems can be 
efficiently approximated by matrix-product states^HMl^ 
in the sense that the value of \ needed to approximate the 
ground state wavefunction to a given accuracy converges 
to a finite value as N — > 00. We therefore think of \ as 
being a finite (but arbitrarily large) number, which can 
be used to control the simulation's precision. 

The matrices T, A can be chosen such tha t they satisfy 
the canonical conditions for an infinite MPS^^Ml 



r A 2 r f =VrU 2 r = i 



(6) 



These equations can be interpreted as stating that the 
transfer matrix 



T aa ';pf3' — 2. r 



/ , x m/3 



AaA 



(7) 



has a right eigenvector Spp/ with eigenvalue A = 1. 
(* denotes complex conjugation.) Similarly, T act i.ppi = 

Sm( r m/3')* r m/3 A " A a' nas a left eigenvector 6 aa t with 
A = 1. We further require that S aa > is the only eigenvec- 
tor with eigenvalue |A| > 1 (which is equivalent to the 
requirement that \tp) is a pure state^). 

The considerations given here become most intuitive 
when one considers, formally, an infinite chain. We form 
a partition of the chain by cutting a bond which results in 
two half-chains. The wavefunction can then be Schmidt 
decomposed^ in the form 



= \ a \aL)\aR), 



(8) 



A. Definitions 

For our simulat ions we use the infinite tree tensor net- 
work state (iTTN)221M29HH] representation of the ground 
state wave function. The iTTN states are the natural 
choice of ansatz state for our model system, since they 
model the tree's geometry. We thus employ this repre- 
sentation to compute the ground state properties numer- 
ically, using the infinite time-evolving block decimation 
(iTEBD) methocP122133| adapted to infinite trees. The 
iTEBD method is a descendant of the densit y matrix 
renormalization group (DMRG) method- ^ 34 * 35 * based on 
matrix product states (MPS^P^l which can be generalized 
to trees. For the sake of completeness, we now review 
some of the properties of MPS's, followed by an intro- 
duction to TTN states. 

A translationally invariant MPS for a chain of length 
L can formally be written in the following form 



\ip) = 2J tr [T mi A. ..T mL A] | mi . ..m L ). 



(5) 



Here, T m are x x X matrices with \ being the dimension of 
the matrices used in the MPS. The index m = —S, . . . , S 



where \aL) and \aR) (a = 1, . . . , x) are orthonormal ba- 
sis vectors of the left and right partition, respectively. In 
the limit of an infinite chain, and under the canonical 
conditions ([Tj]) , the Schmidt values A Q are simply the en- 
tries of the A matrix, A QQ . The A 2 are the eigenvalues 
of the reduced density matrix of either of the two parti- 
tions, and are referred to as the entanglement spectrum. 
The entanglement entropy is S — — J2 a m ^a- This is 
the von Neumann entropy of the reduced density matrix. 
The states \aL) and \aR) can be obtained by multiply- 
ing together all the matrices to the left and right of the 
bond, e.g., if the broken bond is between sites and 

T \*L) = E { , 

Here, 7 is the index of the row of the matrix; when the 
chain is infinitely long, the value of 7 affects only an over- 
all factor in the wavefunction. Reviews of MPS's as well 
as the canonical form can be found in Refs. 1531 and 42. 



B. Tree Tensors Networks 

While MPS are the natural choice of ansatz state for 
one-dimensional systems, they are eminently unsuitable 
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FIG. 1. The two different tilings for the infinite Bethe lattice 
that were employed in the simulations. 

for large higher-dimensional systems. For these systems 
the ideas behind MPS can be generalized to create a 
new class of states known as tensor pr oduct (TPS) or 
projected entangled pair-states (PEPS)P^E31 The con- 
struction of these states is based on the bipartite na- 
ture of trees, allowing them to be split into two subsys- 
tems via the Schmidt decomposition, analogous to one- 
dimensional chains. Thus the generalization of the one- 
dimensional construction to trees is straightforward but 
in order to introduce our ansatz and the notation, we will 
cover it here as well. 

To describe a tree of coordination number z (i.e. each 
vertex has z nearest neighbors), we place tensors rM of 
order z + 1 on the vertices and vectors A fc on the edges 
of the tree graph in Fig. [T] We then connect the tensor's 
indices in a way that mimics the model's underlying lat- 
tice structure. A state on the z = 3 Bethe lattice can 
in this representation be written as 

IVO = (ilfcebonds J2a k < X ^Ofc) X 

(iLgsitcs J2s,<d ^a}a m aSj \ ■ ■ - ) I s i ) I • • ■)• 

While the dimension d of the physical indices is dic- 
tated by the model, the dimension \ of the virtual indices 
afe can be chosen arbitrarily and is only limited by com- 
putational resources. This ansatz can easily be extended 
to lattices with a higher coordination number but in this 
publication wc will only cover the case of z — 3. Analo- 
gous to MPS, the tensors in a tree tensor network can be 
chosen such that they satisfy the conditions for a canon- 
ical tensor network: 

J2a k A2, = 1 (9) 

£ r :Ua m A^A 2 ai ( r:Ua J* = <w (10) 

a k ai 

The advantages of the canonical form of TPS are the 
same as for MPS, i.e. the canonical form provides a 
well-defined basis for evaluations of observables and the 
imaginary-time evolution. 
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FIG. 2. The simple update procedure for a single two-site unit 
cell. By repeating this procedure for every tensor combination 
one full update step is completed. Details on the procedure 
are provided in the text. 



C. Imaginary-time evolution 

In order to obtain the model's ground state within this 
class of iTTN we evolve an initial state \ip) in imagi- 
nary time. Since the Hamiltonian is given by the sum 
over nearest neighbors of products of commuting opera- 
tors, we can implement the imaginary-time evolution as a 
product of local unitary operators using the second-order 
Suzuki- Trotter decomposition of the evolution operator: 
e~ Ht = linin^oo (e~ HSt ) . This decomposition incurs a 
systematic non-accumulating error of 0(St 2 ) which can 
be neglected if the time step St is sufficiently small (in 
our simulations we use 5t = 10 -6 as the final time step). 
While in principle each site tensor and bond tensor A fc 
can be different for different sites and bonds we will cal- 
culate the ground state within the translation-invariant 
sector of iTTN. In this picture, let us consider the ef- 
fects of the imaginary-time evolution of the tree tensor 
network by a translationally invariant Hamiltonian. If 
given a translation-invariant state, the symmetry of the 
Hamiltonian guarantees it to never be broken in time, 
thus enabling us to describe the full state by examining 
only a small number of sites. However, if the system is in 
a phase in which the translational symmetry is sponta- 
neously broken, it turns out to be advantageous to allow 
for a larger unit cell and to guide the wave function into 
a symmetry broken state. 

Since the infinite Bethe lattice is self-similar, a 
translation-invariant ansatz state is the natural choice. 
For numerical reasons it is advantageous to slightly break 
that translational invariance by adopting a larger unit 
cell. Canonically a two-site unit cell is used, but here we 
will extend this scheme and for some calculations employ 
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a six-site unit cell (see Fig. [IJ, which enables us to also 
capture more involved states such as a dimerized state. 
In accordance with the canonical iTEBD algorithm, the 
two-site unit cell consists of two site tensors A and B and 
three bond vectors A x , A v , and A 2 , whereas the six-site 
cell uses six site tensors and nine bond tensors. 

We now have all the ingredients necessary to com- 
pute the ground state. By repeatedly applying the near- 



unitary operators U(St) 



-HSt 



to an initial state \ip) 



and then truncating the entanglement spectrum we can 
obtain the ground state. The operator's near-unitary 
nature allows us to perform the truncation in a well- 
defined basis, yielding a stable algorithm to find the 
ground state on the Bethe lattice. This procedure was 
introduced in the context of MPS as 

iTEBlJiSl 

and later 

on generalized the Bethe latticeP^HH] xhe same algo- 
rithm is also used to find an approximate TPS repre- 
sentation on higher-dimensional lattices in the so-called 
simple update)^ However, as the update algorithm ig- 
nores the loops present in a two-dimensional lattice, the 
TPS found is not optimal to represent the 2D ground 
state. 

The update procedure for the imaginary-time evolu- 
tion of a tree tensor network consisting of an infinitely 
repeated two-site unit cell is now given by the following 
steps (Fig. [2]): 

(1) Contract the site tensors with the adjacent bond 
vectors, leaving one bond open. 

(2) Compute the QR decomposition of the resulting 
tensors relative to the open bond. This modification was 
recently introduced by WangES to reduce the scaling of 
the update with the bond dimension and to stabilize the 
update procedure. 

(3) Contract the evolution operator with the R ten- 
sors resulting from the QR decomposition and calculate 
the singular value decomposition of the resulting 4-index 
tensor. 

(4) Truncate the entanglement spectrum to x entries 
and absorb the unitary matrices in the Q tensors. 

(5) Contract the tensors with the inverse bond vectors 
to obtain the updated site tensors. 

By repeating this procedure for every bond, one full 
update step is completed, i.e., the full lattice is updated. 
Repeatedly applying the update to the system while de- 
creasing the time step brings the initial wave function 
increasingly closer to the true ground state wave func- 
tion. 

For a general Bethe lattice with coordination num- 
ber z the computational cost of this algorithm scales as 
C(x 3z ~ 3 ) + 0(d 2 x z+1 ) and is dominated by the cost of 
the SVD step. The high cost of the SVD limits our al- 
gorithm to values of x ~ 40 on a desktop machine when 
no additional symmetries are used. 



D. Symmetries 

The algorithm's performance can be improved by ex- 
ploiting the model's symmetries, which enable a de- 
composition of the matrices into block-diagonal matri- 
ces, hence reducing the cost of the numerical operations. 
Here, we made use of the models U(l) symmetry in the 
S z sector. By implementing this symmetry the computa- 
tional cost can be significantly reduced, allowing a larger 
cut-off dimension of x ~ 80. Details of how to implement 
this symmetry can be found in e.g. Refs. l47l and l48l 



IV. SYMMETRY PROTECTED TOPOLOGICAL 
PHASES 

Symmetry protected topological phases (SPTP) are 
gapped phases which cannot be characterized by any lo- 
cal order parameters and arc distinct from trivial phases 
(i.e., product states) only in the presence of certain sym- 
metries. In a series of works it had been shown that 
these phases can be completely characterized usi ng pr o- 
jective representations of the symmetries present Eh2li2H31 
In spin systems this means physically that the spin frac- 
tionalizes and the projective representation is due to lo- 
calized spin- half degrees at the edge of a cut. We briefly 
review SPTP for one-dimensional systems and show that 
the concept directly generalizes to the Bethe lattice. 

We start from a state on an infinite chain that is 
invariant under an internal symmetry. The internal sym- 
metry is represented in the spin basis by a unitary matrix 
S acting on each site so that \ip) —> We 
perform a Schmidt decomposition of the system into two 
subsystems (see Fig. |3^) by cutting one bond. We now 
only consider the important Schmidt states which corre- 
spond to Schmidt values A a > e for a given e > 0. 
These Schmidt states transform under a symmetry trans- 
formation as (modulo an overall phase): 



lE(i) 



nR) = ^2u aa ,\a'R), 



(11) 



where U is a unitary matrix which commutes with the A 
matrices Similarly, the left Schmidt states \aL) trans- 
form by the conjugate matrix. As the symmetry element 
g is varied over the whole group, a set of matrices U g 
results. The matrices U g form a x~ dimensional (projec- 
tive) representation of the symmetry group. A projective 
representation is like an ordinary regular representation 
up to phase factors; i.e., if = T, gh , then 



U q U h = e i P^U, 



gh ■ 



(12) 



The phases p{g, h) can be used to classify different topo- 
logical phasesPEHU 

As a concrete example, we now consider the Hal- 
dane phase around the AKLT state in the presence of 
a Z 2 x Z 2 symmetry. The generators of the symme- 
try group are the spin rotations 1Z X — exp(iitS x ) and 
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= 7i(in>-ut» 



(b) 




singlet and thus the Schmidt states have localized half- 
integer spins at the edges. Thus the transformation of 
the Schmidt states under a symmetry operation yields 
a (projective) representation U of the symmetry group 
which characterizes the phase. Note that there are in 
fact different kinds of SPTP that can be realized on the 
Bethe lattice depending on the fractionalization of the 
spin. For example, in an S = 1 model on a z = 3 lat- 
tice with strong dimerization (two strong bonds and one 
weak bond) , we would obtain a network of S = 1 Haldane 
chains. 




FIG. 3. Diagrammatic representation of the AKLT wave 
function on a chain (a) and a z — 3 Bethe lattice (b). The 
red ovals are representing S — 1/2 singlets. A Schmidt de- 
composition at the blue line cuts in both cases one singlet, 
leaving behind localized S = 1/2 edge spins. 



1Z Z = exp(inS z ). The phases for each spin rotation indi- 
vidually (e.g., U% = e ta t) can be removed by redefining 
the phase of the corresponding [/-matrix. However, the 
representations of 1Z X 1Z Z and 1Z Z 1Z X can also differ by a 
phase, which it turns out must be ±1: 



u x u z = ±u z u x 



(13) 



I.e., the matrices either commute or anti-commute. This 
resulting phase cannot be gauged away because the 
phases of U x and U z enter both sides of the equation 
in the same way. Thus we have two different classes 
of projective representations. If the phase is —1, then 
the spectrum of A is doubly degenerate, since A com- 
mutes with the two unitary matrices U x , U z which anti- 
commute among themselves. For the AKLT state consid- 
ered here, the Schmidt states have half-integer edge spins 
(see Fig. [3]). Thus we find U x — a x and U z = a z , there- 
fore U X U Z = —U Z U X , and the Haldane phase is protected 
if the system is symmetric under both 1Z X and 1Z Z . An 
analogous argument can be made for inversion symmetry 
(i.e., spatial inversion of the system at a bond) and time 
reversal symmetry 

Using the above arguments, we can now characterize 
the 5 = 3/2 AKLT state on the z = 3 lattice. As il- 
lustrated in Fig. [3)3, a single cut through a bond sepa- 
rates the Bethe lattice into two disconnected subsystems. 
In the AKLT state, any of the bonds has a S = 1/2 



V. NUMERICAL RESULTS 

We employ our variant of the iTEBD algorithm to 
study the bilinear-biquadratic-bicubic Heisenberg model 
defined in Eq. ([!]) over the full parameter range of 
ip = —ix . . . 7r, 8 = To ensure unbiased re- 

sults we use different iTEBD implementations, one with 
a two-site unit cell and one with a six-site unit cell (see 
Fig. [lj. The simulations were conducted both with and 
without explicitly conserving the model's U(l) symme- 
try and were started from different initial wavefunctions 
ranging from completely random to fully polarized initial 
states. Furthermore, we also studied the dependence on 
the evolution scheme by starting the imaginary-time evo- 
lution with a slightly modified evolution operator (e.g. by 
adding a small ferromagnetic or symmetry-breaking in- 
teraction term to the Hamiltonian) and only later in the 
calculation using the actual Hamiltonian to converge the 
trial wavefunction to the groundstate. 

The simulations were started with a large time step 
of St = 0.1 which was then in several steps decreased 
to 5t — 1CP 6 , reducing the Trotter error to insignifi- 
cance. These checks are commonly accepted best prac- 
tice for any imaginary-time evolution, as the overlap of 
the initial state and the true groundstate has to be fi- 
nite in order for the algorithm to be able to evolve the 
trial wavefunction to the groundstate. Hence we have 
to ascertain the result's independence of the procedure 
and the initial states. As expected we find the strongest 
dependence on the initial tensors in the ferromagnetic 
(anti-ferromagnetic) phase when starting with an anti- 
ferromagnetic (ferromagnetic) initial state, but no signif- 
icant correlation otherwise. To establish that our findings 
are not dependent on the unit cell, i.e., the constraint of 
a two-site unit cell, we also implemented a variant of 
the iTEBD algorithm that operates on a six-site unit 
cell, see Fig. 1(a) This implementation enables us to 



also describe a dimerized phase which would otherwise 
be able to escape characterization as it is not commen- 
surable with a two-site unit cell. 

Special attention was paid to the four S£/(4) symmet- 
ric points shown in the phase diagram Fig. [4] where 
we ran a variety of simulations to ascertain the model's 
behavior. With the exception of the multi-critical point 
at <p = 2.93, 8 = 0.17 we were not able to observe any 
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FIG. 4. The phase diagram as determined by the tree net- 
work iTEBD method. We established the presence of an anti- 
ferromagnetic phase (7), a ferromagnetic phase (77), the Hal- 
dane phase (777) , a dimerized phase (IV) , and a ferrimagnetic 
phase (V). For reference purposes the phase diagram also in- 
cludes the AKLT point denoted by a black star and the SO (5) 
symmetric line denoted by the dotted line between phases 7 
and 77. The SU(4) symmetric points are marked by white 
stars and lie on the SO (5) line. 



FIG. 5. The (staggered) magnetization for various of values 
of and ip — 0.9 with \ — 24. The ferromagnetic phase 
can be identified as the region with maximal uniform mag- 
netization per site, whereas the anti-ferromagnetic phase is 
distinguished by the finite staggered magnetization and van- 
ishing uniform magnetization per site. Inbetween those two 
phases are the dimerized and the Haldane phases, where both 
the staggered (mj) and uniform magnetization (m z ) vanish. 
In the ferrimagnetic phase they take intermediate values. The 
colors correspond to the phases introduced in Fig. [4] 



behavior diverging from the encompassing phase. 

Our results for the phase diagram of the bilinear- 
biquadratic-bicubic Heisenberg model are shown in 
Fig. [4] We find five different phases: a Neel phase, a fer- 
romagnetic polarized phase, the Haldane phase, a dimer- 
ized phase, and a ferrimagnetic phase. 



A. (Anti-)Ferromagnetic phases 

We start our analysis of the phase diagram with the 
straightforward phases, i.e. the AFM and FM phases. 
First we consider the staggered magnetization m s z and 
the uniform magnetization m z . The polarized phases 
can easily be identified by observing the two different 
order parameter's behavior shown in Fig. [5] In the fer- 
romagnetic region 8 < — | the staggered magnetization 
m s z disappears, while the magnetization per site m z is 
maximal in this region. The opposite holds for the anti- 
ferromagnetic part of the phase diagram, where only the 
staggered magnetization remains non-zero. Inspired by 
the (anti-) nematic order found in mean-field studies^ 
we also calculate the octupolar order parameter in the 
vicinity of the SO (5) line. Again we observe the absence 
of different order, except for some small contributions in 
the ferromagnetic region at negative ip. This might be 
due to remnants of the classical order persisting to zero 
temperature. 
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FIG. 6. Behavior of the dimer order parameters defined in 
Eq. [14] for a cut through the phase diagram at ip = 0.9 with 
X = 24. The increase of the dimer order parameters occurs 
at the boundary between the Haldane and the ferrimagnetic 
phase. 



B. Dimerized phase 

We determine a dimerized phase to be present in the 
region denoted by IV in the phase diagram (Fig. Q . This 
phase is hard to characterize due to its diminutive size 
and vanishing magnetization, however, careful calcula- 
tions strongly indicate its existence. To determine the 
properties of this elusive plains we calculate the xy and 
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FIG. 7. The entropy for a cut in 8 direction through the 
phase diagram at tp = 0.9 for various values of \ showing 
the entropy to be clearly converged when increasing the bond 
dimension. 



z components of the dimer order parameters D xy and 
D xy defined via 

iK'U ((sfsy + sfs*) (i4) 

-(Sf + (15) 

DU k = (s?s;-sjsd, (i6) 

where i, j, and k label consecutive lattice sites residing 
on different shells (e.g. sites A, B, and C in Fig. [l](b)). 
Calculation of these order parameter components reveals 
them to vanish for the magnetically ordered phases, as 
well as in the Haldane phase. Only in the dimer phase 
do they assume finite values. The hypothesis of a dimer 
phase is further corroborated by the vanishing magne- 
tization. As opposed to the Haldane phase we also fail 
to observe finite edge spins in this phase. Together these 
observations indicate the existence of a narrow dimcrized 
phase close to the 50(5) symmetric line. 



C. Ferrimagnetic phase 

We found a ferrimagnetic phase that exists between the 
dimer phase and the ferromagnetic phase. This phase 
shows both a finite staggered magnetization as well as 
a finite magnetization per site, but displays vanishing 
dimer order parameters. As a test for this phase we try 
adding a small (staggered) field in the z direction to the 
Hamiltonian. Performing the simulation with this mod- 
ified Hamiltonian results in strong polarization which 
could be both anti-ferromagnetic or ferromagnetic de- 
pending on the applied field. 
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FIG. 8. Entanglement spectrum throughout a cut in the 
phase diagram at cp — 0.9 with x — 24 for various values of 
6. Clearly visible is the doubling of the levels in the Haldane 
phase. 



D. Haldane phase 

The last phase is a phase with exponentially decaying 
correlations and no broken symmetries (i.e., no local or- 
der parameter exists) which is identified as a SPTP. In 
this phase, the entanglement spectrum always displays 
an even degeneracy in the entire phase, which is clearly 
visible in the iTEBD calculations. The identification of 
this phase rests on the presence of a finite 5 = 1/2 edge 
spin and its characteristic degeneracies in the entangle- 
ment spectrum shown in Fig. [HJ As a direct evidence, 
we also calculated the edge spin of the Schmidt states 
directly and find a localized spin-1/2. In our simula- 
tions we are also able to observe the AKLT point where 
the ansatz state reduces to an exact tensor network with 
bond dimension \ — 2- By adding small perturbations 
which destroy all the necessary symmetries to protect 
the phase, it is possible to drive the system out of the 
Haldane phase without a phase transition (i.e., the de- 
generacies in the spectrum are lifted for an arbitrarily 
small perturbation). We also checked numerically that 
the phase is robust against small perturbations which do 
not break the symmetries needed to protect it. All of 
these observation can be explained by the presence of a 
SPT phase as discussed in section [TV] 

VI. SUMMARY 

In this work, we have studied a general SU(2) sym- 
metric spin-3/2 model on the z = 3 Bcthe lattice using 
the infinite Time Evolving Block Decimation (iTEBD) 
method. We found that the model exhibits a rich phase 
diagram containing several magnetic phases, a dimer- 
ized phase as well as a symmetry protected topologi- 
cal phase (SPTP). The magnetic phases were identified 



by calculating the uniform and the staggered magneti- 
zation. We found a polarized ferromagnetic phase, an 
anti-ferromagnetic as well as a ferrimagnetic phase with 
finite (staggered) magnetization. Our simulations sug- 
gest the presence of a dimerized phase with vanishing 
magnetization and finite dimer order parameters. We 
also identified a symmetry protected topological phase 
which shows all the key features of the Haldane phase. 
This phase is characterized by spin- 1/2 edge spins and 
degeneracies in the entanglement spectrum. 
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